using System;
using System.Collections;
using System.Collections.Generic;
using NUnit.Framework;
using dnAnalytics.Math;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using dnAnalytics.LinearAlgebra.IO;
using dnAnalytics.Exceptions;

namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
{
    public abstract class AbstractLUTest
    {
        public abstract void GetDecompositionClass(String file, out Matrix matrix, out LU lu);
        private LU lu;
        private LU singularLu;
        private Matrix matrix;
        private Matrix singularMatrix;

        public abstract Matrix GetMatrix(int order);
        public abstract Matrix GetMatrix(int rows, int columns);
        public abstract Vector GetVector(int size);

        [TestFixtureSetUp]
        public void SetUp()
        {
            MatrixMarketReader mmr = new MatrixMarketReader();
            GetDecompositionClass("../../TestMatrices/random_real_general_array_100.matrix_market", out matrix, out lu);
            GetDecompositionClass("../../TestMatrices/gear_integer_general_coordinate_100.matrix_market", out singularMatrix, out singularLu);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void LUConstructorNull()
        {
            Matrix matrix = null;
            LU lu = new LU(matrix);
        }

        [Test]
        [ExpectedException(typeof(NotSquareMatrixException))]
        public void LUConstructorNotSquare()
        {
            Matrix matrix = GetMatrix(3, 2);
            LU lu = new LU(matrix);
        }

        [Test]
        public void Singular()
        {
            Assert.IsTrue(singularLu.IsSingular());
        }

        [Test]
        public void NotSingular()
        {
            Assert.IsFalse(lu.IsSingular());
        }

        [Test]
        public void Factors()
        {
            int[] pivots = lu.Pivots();

            Matrix perm = GetMatrix(pivots.Length);
            for (int i = 0; i < pivots.Length; i++)
            {
                perm[i, i] = 1;
            }

            //there has to be a quicker way to generate the perm squareMatrix?
            double temp;
            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    for (int j = 0; j < pivots.Length; j++)
                    {
                        temp = perm[i, j];
                        perm[i, j] = perm[pivots[i], j];
                        perm[pivots[i], j] = temp;
                    }
                }
            }

            Matrix lower = lu.LowerFactor();
            Matrix upper = lu.UpperFactor();
            Matrix result = GetMatrix(lower.Rows);
            lower.Multiply(upper, result);

            Matrix expected = GetMatrix(lower.Rows);
            perm.Multiply(matrix, expected);

            double error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(expected[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - expected[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - expected[i, j]) / expected[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }


            pivots = singularLu.Pivots();
            perm = GetMatrix(pivots.Length);
            for (int i = 0; i < pivots.Length; i++)
            {
                perm[i, i] = 1;
            }

            //there has to be a quicker way to generate the perm squareMatrix?
            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    for (int j = 0; j < pivots.Length; j++)
                    {
                        temp = perm[i, j];
                        perm[i, j] = perm[pivots[i], j];
                        perm[pivots[i], j] = temp;
                    }
                }
            }

            lower = singularLu.LowerFactor();
            upper = singularLu.UpperFactor();
            result = GetMatrix(lower.Rows);
            lower.Multiply(upper, result);

            expected = GetMatrix(lower.Rows);
            perm.Multiply(singularMatrix, expected);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(expected[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - expected[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - expected[i, j]) / expected[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        public void FactorsResults()
        {
            int[] pivots = lu.Pivots();
            Matrix perm = GetMatrix(pivots.Length);
            for (int i = 0; i < pivots.Length; i++)
            {
                perm[i, i] = 1;
            }

            //there has to be a quicker way to generate the perm squareMatrix?
            double temp;
            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    for (int j = 0; j < pivots.Length; j++)
                    {
                        temp = perm[i, j];
                        perm[i, j] = perm[pivots[i], j];
                        perm[pivots[i], j] = temp;
                    }
                }
            }

            Matrix lower = GetMatrix(matrix.Rows);
            Matrix upper = GetMatrix(matrix.Rows);

            lu.LowerFactor(lower);
            lu.UpperFactor(upper);
            Matrix result = GetMatrix(lower.Rows);
            lower.Multiply(upper, result);

            Matrix expected = GetMatrix(lower.Rows);
            perm.Multiply(matrix, expected);

            double error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(expected[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - expected[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - expected[i, j]) / expected[i, j]);
                    }

                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }


            pivots = singularLu.Pivots();
            perm = GetMatrix(pivots.Length);
            for (int i = 0; i < pivots.Length; i++)
            {
                perm[i, i] = 1;
            }

            //there has to be a quicker way to generate the perm squareMatrix?
            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    for (int j = 0; j < pivots.Length; j++)
                    {
                        temp = perm[i, j];
                        perm[i, j] = perm[pivots[i], j];
                        perm[pivots[i], j] = temp;
                    }
                }
            }

            lower = GetMatrix(matrix.Rows);
            upper = GetMatrix(matrix.Rows);

            singularLu.LowerFactor(lower);
            singularLu.UpperFactor(upper);
            result = GetMatrix(lower.Rows);
            lower.Multiply(upper, result);

            expected = GetMatrix(lower.Rows);
            perm.Multiply(singularMatrix, expected);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(expected[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - expected[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - expected[i, j]) / expected[i, j]);
                    }

                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void LowerFactorResultNull()
        {
            Matrix result = null;
            lu.LowerFactor(result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void UpperFactorResultNull()
        {
            Matrix result = null;
            lu.UpperFactor(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LowerFactorResultNotConformable1()
        {
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Rows);
            lu.LowerFactor(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void LowerFactorResultNotConformable2()
        {
            Matrix result = GetMatrix(matrix.Rows, matrix.Rows + 1);
            lu.LowerFactor(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void UpperFactorResultNotConformable1()
        {
            Matrix result = GetMatrix(matrix.Rows + 1, matrix.Rows);
            lu.UpperFactor(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void UpperFactorResultNotConformable2()
        {
            Matrix result = GetMatrix(matrix.Rows, matrix.Rows + 1);
            lu.UpperFactor(result);
        }

        [Test]
        public void Determinant()
        {
            Assert.AreEqual(0, singularLu.Determinant());
            double expected = -5.006656015927541e248;
            double actual = lu.Determinant();
            double error = System.Math.Abs((actual - expected) / expected);
            Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
        }

        [Test]
        public void Inverse()
        {
            Matrix inverse = lu.Inverse();
            Matrix result = GetMatrix(inverse.Rows);
            matrix.Multiply(inverse, result);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void InverseSingular()
        {
            Matrix inverse = singularLu.Inverse();
        }


        [Test]
        public void InverseResult()
        {
            Matrix inverse = GetMatrix(matrix.Rows);
            lu.Inverse(inverse);
            Matrix result = GetMatrix(inverse.Rows);
            matrix.Multiply(inverse, result);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void InverseResultSingular()
        {
            Matrix inverse = GetMatrix(matrix.Rows);
            singularLu.Inverse(inverse);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void InverseResultNull()
        {
            Matrix inverse = null;
            lu.Inverse(inverse);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void InverseResultNotConformable1()
        {
            Matrix inverse = GetMatrix(matrix.Rows + 1, matrix.Rows);
            lu.Inverse(inverse);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void InverseResultNotConformable2()
        {
            Matrix inverse = GetMatrix(matrix.Rows, matrix.Rows + 1);
            lu.Inverse(inverse);
        }

        [Test]
        public void L1NormConditionNumber()
        {
            double expected = 0.000539909549815265;
            double actual = lu.L1NormConditionNumber();
            double error = System.Math.Abs(actual - expected);
            Assert.IsTrue(error < Constants.ACCEPTABLE_L1NORM_CONDITION_NUMBER_ERROR);
        }

        [Test]
        public void InfinityNormConditionNumber()
        {
            double expected = 0.000522256126360585;
            double actual = lu.InfinityNormConditionNumber();
            double error = System.Math.Abs((actual - expected) / expected);
            Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
        }

        [Test]
        public void SolveMatrix()
        {
            Matrix identity = MatrixBuilder.CreateIdentityMatrix(matrix.Rows);
            Matrix result = lu.Solve(identity);
            result = matrix * result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveMatrixSingularMatrix()
        {
            Matrix identity = GetMatrix(matrix.Rows);
            singularLu.Solve(identity);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNull()
        {
            Matrix identity = null;
            lu.Solve(identity);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformable()
        {
            Matrix identity = GetMatrix(matrix.Rows + 1);
            lu.Solve(identity);
        }

        [Test]
        public void SolveMatrixResult()
        {
            Matrix identity = MatrixBuilder.CreateIdentityMatrix(matrix.Rows);
            Matrix result = GetMatrix(matrix.Rows);
            lu.Solve(identity, result);
            result = matrix * result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveMatrixSingularMatrixResult()
        {
            Matrix identity = GetMatrix(matrix.Rows);
            Matrix result = GetMatrix(matrix.Rows);
            singularLu.Solve(identity, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult1()
        {
            Matrix identity = null;
            Matrix result = GetMatrix(matrix.Rows);
            lu.Solve(identity, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult2()
        {
            Matrix identity = GetMatrix(matrix.Rows);
            Matrix result = null;
            lu.Solve(identity, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult1()
        {
            Matrix identity = GetMatrix(matrix.Rows + 1);
            Matrix result = GetMatrix(matrix.Rows);
            lu.Solve(identity, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult2()
        {
            Matrix identity = GetMatrix(matrix.Rows);
            Matrix result = GetMatrix(matrix.Rows + 1);
            lu.Solve(identity, result);
        }

        [Test]
        public void SolveVector()
        {
            Vector vector = VectorBuilder.CreateVector(matrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (double)(i + 1) / (i + 2);
            }
            Vector result = lu.Solve(vector);
            Vector resultVector = matrix * result;
            for (int i = 0; i < resultVector.Count; i++)
            {
                Assert.AreEqual((double)(i + 1) / (i + 2), resultVector[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveVectorSingular()
        {
            Vector vector = GetVector(matrix.Rows);
            singularLu.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNull()
        {
            Matrix identity = null;
            lu.Solve(identity);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformable()
        {
            Matrix identity = GetMatrix(matrix.Rows + 1);
            lu.Solve(identity);
        }

        [Test]
        public void SolveVectorResult()
        {
            Vector vector = VectorBuilder.CreateVector(matrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (double)(i + 1) / (i + 2);
            }

            Vector result = VectorBuilder.CreateVector(matrix.Rows);
            lu.Solve(vector, result);
            Vector resultVector = matrix * result;
            for (int i = 0; i < resultVector.Count; i++)
            {
                Assert.AreEqual((double)(i + 1) / (i + 2), resultVector[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        [ExpectedException(typeof(SingularMatrixException))]
        public void SolveVectorSingularResult()
        {
            Vector vector = GetVector(matrix.Rows);
            Vector result = GetVector(matrix.Rows);
            singularLu.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult1()
        {
            Vector vector = null;
            Vector result = GetVector(matrix.Rows);
            lu.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult2()
        {
            Vector vector = GetVector(matrix.Rows);
            Vector result = null;
            lu.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResult1()
        {
            Vector vector = GetVector(matrix.Rows + 1);
            Vector result = GetVector(matrix.Rows);
            lu.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResult2()
        {
            Vector vector = GetVector(matrix.Rows);
            Vector result = GetVector(matrix.Rows + 1);
            lu.Solve(vector, result);
        }
    }
}
